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V 

surface variable along orthogo- 
nal trajectory of rays 

e 

angle, relative to N , at which 
ray is launched from surface 

K g i 

geodesic curvature of ray 

k 9 2 

geodesic curvature of orthogo- 
nal trajectory of rays 

Kl 

component of curvature 
vector of surface wave front 
in direction p 

K n 1 

normal curvature of ray 

K n 2 

normal curvature of orthogonal 
trajectory of rays 

Kyj 

curvature vector of surface 
wave front 

X 

direction ^ on surface 

A n 

direction ^ orthogonal to rays 

P> 

general ray parameter 

£ 

specific surface parameter 

P 

parameter determined by 
densities of materials on two 
sides of surface 

a 

arc length along ray 

<t> 

II 

1 

$ 

angle that surface paraxial ray 
makes with N,T plane 

Subscripts: 


/ 

base point of surface 

i 

initial point of surface 

m 

conditions in medium sur- 
rounding surface 

0 

initial point on surface ray 

Basic Considerations 

The analysis utilizes vector terminology consis- 
tent with such fundamental differential geometry 
texts as references 7 and 8. The section “Analysis” 
includes a derivation of the required equations and a 
description of the procedure for computing the sur- 
face rays and the ray strip width, which is analogous 
to ray tube area for volume waves. These quantities 


are required to obtain the phase and amplitude of 
the surface wave. 

For purposes of analysis, it is assumed that the 
surface can be described by a vector equation in two 
surface variables u, v. That is, a point on the surface 
is given by 

r(u, u) = x(u, v) i + y(u , if) } T z(u, v) k (1) 

A class of practical shapes that can be modeled 
in this form is discussed in the section “Analytically 
Lofted Surfaces.” The section “Sample Calculations” 
illustrates how the wave equations are used to com- 
pute surface waves on such surfaces. The final section 
treats some details of computing the wave field and 
considers the problem of waves launched from the 
surface into a surrounding medium that will support 
waves. 

The analysis does not include derivation of the 
wave speeds associated with the various types of 
surface waves. These derivations are included in 
the basic works on those waves, as for example, 
reference 9. If the surface is such that the wave speed 
is uniform, it follows from Fermat’s principle that the 
surface wave propagates along geodesic lines; that is, 
the ray paths are geodesics. 

If the wave speed is not uniform but is given 
as a function of position, the rays no longer follow 
geodesic paths; however, the surface wave fronts are 
still orthogonal trajectories of the rays. Problems 
for which the wave speed is determined as a function 
of direction are not treated. Important problems in 
this category include surface waves on anisotropic 
crystals, Rayleigh waves on surfaces having local 
radii of curvature not much larger than the wave 
length, and certain Lamb waves in relatively thick 
shells with relatively small radii of curvature. 

Analysis 

With the surface described in terms of the surface 
variables u,v by equation (1), the element of arc 
length is given by 

ds 2 = dr • dr = E dv? + 2 F du dv + G dv 1 (2) 
where E , F, G are the metric coefficients 

E — v u * v u (3a) 

F = r u • r v (3b) 

G = v u ■ v v (3 C ) 

A direction on the surface is determined by assigning 
a value A to the derivative ^ . Thus, if an arbitrary 
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increment du is taken, the corresponding arc incre- 
ment ds in the A direction is 


ds = 


\/ £+2f i +G W du 

\/e + 2FX + GA 2 du 


(4) 


The u, v variables are the basic geometric param- 
eters used to define the surface, and calculations are 
normally performed in this system. However, since 
the wave fronts are orthogonal trajectories of the 
rays, it is convenient analytically to utilize an alter- 
nate set of coordinates taken, respectively, along the 
rays and along the orthogonal trajectories. Along the 
rays, we take time t as the coordinate and denote arc 
length along the ray by a so that 


d<7 — c dt — \fW 0 dt (5) 


The coordinate normal to the rays is denoted by 
77 . Thus 77 determines which ray of the family is 
being considered. (See fig. 1 .) Denote by dn the 
incremental arc length normal to the rays: 

dn = \J~Go dv ( 6 ) 


the direction normal to ray. This is an intrinsic rela- 
tion and therefore is true in any coordinate system. 
For purposes of calculation it must be expressed in 
terms of the u, v coordinates that define the surface. 
Toward this end, assume first that an initial point 
u 0 ,v Q is known on the ray, along with an initial di- 
rection X 0 and the distribution of wave speed c(u, v). 

The direction normal to the ray is obtained from 
the relation 



(See ref. 7, p. 59, eq. (2-10).) Then ^ is obtained 
from equation (4): 


du _ 1 

dn + 2FX n + G\l 


( 11 ) 


and 

dv _ 6v du __ ^ du 
dn 6udn n dn 

With equations (11) and (12), the derivative 


( 12 ) 


dc _ dc du dc dv 
dn du dn dv dn 


(13) 


Thus, the general arc length element is 

ds 2 = E 0 dt 2 + G 0 dt] 2 (7) 

where E 0 = c 2 . 

We could obtain a differential equation for the 
rays by minimizing the time integral 

t-t\= I — ( 8 ) 

J(T\ C 

utilizing a procedure analogous to that used to de- 
rive the equations of a geodesic by minimizing the 
distance integral (ref. 7, p. 140). However, a briefer 
and simpler approach is available because of the or- 
thogonality of the £, 77 coordinate system. It is known 
(ref. 7, p. 130) that the geodesic curvature of the first 
set of coordinate lines (77 = Constant) is given by the 
formula 

1 dE o 

Kgl ~ 2E 0 \fU~ 0 drj 


This equation states that the geodesic curvature 
of a ray is proportional to the derivative of log c in 


in equation (9) is expressed in terms of u and v. 
The geodesic curvature K g \ can also be expressed in 
terms of the original coordinates (u,u) and metric 
coefficients. This relation is (see ref. 7, p. 128, 
eq. ( 1 - 6 )) 

= [In(t0 3 + ( 2 r? 2 - rioMV + (rh - 2 ri 2 K(t/) 2 


— r \ 2 {v') i + (u'v"-u"v')\H (14) 

where primes indicate differentiation with respect to 
a, and 

H = s/EG - F 2 (15) 

The T’s are Christoffel symbols defined as follows 
(ref. 7, p. 107) 

2tf 2 r} 1 = GE U - 2 FF U + FE V (16a) 

2tf 2 r} 2 = GE V - FG U ( 16 b) 

2tf 2 r^ = 2 GF V - GG U - FG V (16c) 

2H 2 Y 2 n = 2EF U - EE V - FE U (16d) 

2H 2 T\ 2 = EG U - FE V (16e) 

2H 2 r% 2 = EG v -2FF v + FG u (16f) 
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The differential equation for the ray path in the 
u,v system is now obtained by substituting rela- 
tion (14) into equation (9). To simplify the resulting 
equation, we first change the ray parameter from arc 
length a to a general ray parameter p, This is ac- 
complished by multiplying through by (jfy : 



+ (Yl 2 -2T\ 2 )u\v'j 1 


- r£ 2 (t/) 3 + ( uV' 




1 dc 
c dn 


(17) 


Now, since 

f dc r\ 3 du d?v __ du d 2 v 
\dp ) dcr da 2 dp dp? 

and similarly 


f da\^ d?u dv _ d?u dv 
\dp J da 2 da dp 2 dp 


(18a) 


(18b) 


it follows that equation (17) can be written as 


r 2 n («') 3 + ( 2r } 2 - r },)( u')V + (t 2 22 - 2 t\ 2 )u'(v ') 2 - r » 2 (*/) 3 


-I- (u u — u'v) 


»"(£)*;£ < 19 > 


where now the primes denote differentiation with 
respect to p. Finally, we choose for the arc 
parameter the coordinate variable u, so that 
u f = \,u n — 0, v 9 = = A, the local ray direc- 

tion. With these substitutions, equation (19) can be 
written as an equation for the local change in ray 
direction. 


increment in v: dv = A du, iq = v 0 + dv; and 
then computing c, and the T’s at u\,v\. 

Substituting these quantities into the right-hand side 
of equation (20) and multiplying by du yields the 
change in direction dX for the next step. Then A 
is incremented by this amount and a step is taken 
in the new direction; this process is continued to 
generate the ray. The arc length on each step is 
calculated from the elapsed time (eq. (5)). Accuracy 
of the procedure can be checked by duplicating a 
calculation with different step sizes. Computation 
time for 1000 steps is at most a few seconds on a 
modern computer. 

To determine the amplitude variation, the spread- 
ing of the wave must be calculated. This spreading 
is measured by the distance between the rays: dn , in 
equation (6). Since for a pair of rays, drj is fixed and 
since dn = \[G~ 0 drj, the problem is to determine G 0 
as a function of a. 

When the surface does not have a uniform wave 
speed, the simplest way to approach this problem, 
both conceptually and analytically, is to perform the 
calculation directly. We simply compute two adja- 
cent rays simultaneously. Equal increments in time 
are specified, so that in equation (20) the increments 
du are determined by the relation 


du(t) = 


du da 
d^~dt dt 


c dt 

s/E + 2FA + GA 2 


( 22 ) 


Thus, at each step we have a point u\,v\ on the first 
ray and a corresponding point u<i,V 2 on the second 
so that 

dn(t) = ^(Au) 2 + 2F Au Ar; +G( Au) 2 (23) 


^ — r£ 2 A 3 + (2r} 2 - r 22 )A 2 + (r}j - 2r 2 2 )A - 


1 dc f da \ 3 
cH dn \du ) 


( 20 ) 


where, by equation (4), 


^ = s/E + 2FA + GX 2 (21) 

du 

Equation (20) is used to compute a ray as follows. 
Starting at the initial point u Q , v Q one takes a step in 
the initially specified direction X Q by incrementing 
u 0 : u\ = u Q + du; computing the corresponding 


and 

\/ g » (,) - W <24) 

where dr] is the constant increment distinguishing the 
two rays. 

It is possible to compute the ray strip width 
without calculating an adjacent ray independently. 
Since this kind of development may be required for 
analytical purposes, the theory will be outlined here 
as a matter of academic interest. 

The variation of y/G^ with t , is related to 

L and to the local total curvature K by Gauss’s 
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equation, which for this case, is (ref. 7, p. 113) 


d 2 VG^ 1 d?E 0 . 

w~ + w~^* ^ (25) 


where the total, or Gaussian, curvature K is 


geodesic curvature n g 2 of the wave fronts by an equa- 
tion similar to equation (14). The ray spreading %/C?^ 
is related to K g 2 by the formula 

_ 1 d\/Gl _ 1 dy/G^, /01 ^ 

Kg2 ~ y/TT 0 da “ Cs [Gl dt 1 } 


K = 


eg - / 2 
H 2 


(26) 


The quantities e, /, g are the curvature coefficients or 
coefficients of the second fundamental form. They 
are computed by the formulas 


— Tun * N 1 


/ = ~~ ’ N 

%uv 


9 = — • N 

X V v 


(27) 


where N is the local surface normal 


N = 


r u x r„ 

H 


(28) 


Thus, K can be computed at each point of the ray, 
and since \[Eo = c, equation (25) can be written 


^ + (29) 

For a surface with uniform wave speed, the right- 
hand side of equation (29) vanishes, and \fG~ 0 can 
be calculated in a straightforward step-by-step man- 
ner as the points on the ray (a geodesic in this 
case) are computed. (Application of the WKB 
method (ref. 10) would involve complications, since 
K changes sign when the ray traverses an inflected 
region.) In this case, this procedure may well be 
preferable to the alternate method of computing an 
adjacent ray. 

However, if the right-hand side of equation (29) 
does not vanish, we have 


d 2 c _ d 2 c ( du \ 2 d 2 c /dt >\ 2 dc d 2 c dc d 2 v 

dn 2 du 2 \dnj + dv 2 ydny + du dn 2 + dv dn 2 

(30) 

The derivatives of c with respect to u and v can 
be calculated from the given distribution e(u, v). 
The first derivatives ^ are calculated by equa- 
tions (10) through (12). The problem lies with the 

second derivatives, They are related to the 


(ref. 7, p. 130, eqs. (1)-(10)). 

Thus, given the local geodesic curvature of the 
wave front at the initial point, one can compute 

^ and in a manner similar to that applied to 
equation (14). Then equation (29) can be integrated 

one step to obtain and \[G~ 0 at the next point. 
Next, K g 2 can be calculated by equation (31), and 
thus the process is continued. Thus, it is seen that, 
with a nonuniform wave speed, this procedure is more 
difficult than simply computing an adjacent ray in 
theory, in programming, and in calculation time. 

Analytically Lofted Surfaces 

The procedure described in the previous section 
can be applied to any surface that is given in the form 
of equation (1). The example given here will utilize 
a class of analytically lofted surfaces described in 
reference 6. These surfaces are defined by specifying 
the shape of an initial cross section in a plane x = 
and a final, or terminating, cross section shape at x — 
x f and then requiring the intermediate cross section 
to represent a smooth transition from the initial to 
the final shape. The shapes of the two extreme 
cross sections y(£), z(£) are described parametrically 
in terms of an arbitrary parameter £. This parameter 
might be for example the normalized y coordinate, 
the normalized arc length, the angle variable in polar 
coordinates, or the ellipse angle parameter. 

Let a(x) be a function, varying smoothly and 
monotonically from 0 to 1 as x varies from X{ to xj 

(homotopy function). Examples are and 

HiSi)]" Then, for Xi < x < xj a cross- 
section shape is defined by 

y{x, 0 = [1 - a(*)] Vi(0 + a(x) yj{t) (32a) 
z(x, 0 = [1 - a(®)] Zi(£) + a(x) Zf(£) (32b) 

These shapes vary smoothly from the initial shape 
at X{ to the terminating shape at xj. By taking 
£ as a normalized coordinate, the size variation of 
the cross sections can be specified by prescribing a 
scaling function /?(#), which varies smoothly but not 
necessarily monotonically from 0 to 1 as x varies from 
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X{ to xj. The scaled variables are defined by 


the lofting functions were taken as 


v(x,0 = P(x) y(x,Z) (33a) 

z{x, 0 = P(x) z(x,Z) (33b) 

The vector equation of this transition surface is 

r{x,0 = xi + y(x,Oj + z(x,Ok ( 34 ) 

which has the required form (eq. (1)) where the sur- 
face variables u,v are x,£, respectively. For model- 
ing complicated shapes, such as blended wing-body 
combinations, 0 must be allowed to vary with f as 
well as x, but for illustrative purposes equations (33) 
describe a sufficiently general class of surfaces. The 
x = Constant lines describe the cross-section shapes, 
and the £ = Constant lines are the lofting lines. 

The first derivatives are 


d/3 ^ da 
r ’ =,+ \di v+0 di 
d/3 . da 


■*i) 


(35a) 


r$ = 0 

+ 0 


n , d Vi dy f 
dzi dzj 




(35b) 


From these formulas, the metric coefficients can be 
calculated as well as H and, if needed, the surface 
normal N. Equations (35) can be differentiated again 
if the curvature coefficients are required. 


Sample Calculations 

For purposes of illustration, a surface is defined 
as follows. The initial and final cross sections are 
defined to be ellipses: 


With 


Vi(0 = cos£ 
Zi(0 = ^sinf 
?//(£) = cos£ 
z f{0 = 0.58 sin £ 


R{x) 


X-^Xj 

Xf X i 


a{x) = R 

«*) = r 3/4 (i - fi) 1/2 

A top view of this surface is shown in figure 2. 

To illustrate the wave calculation, we assume a 
source located on the surface at 


x 0 = 0.55 (xj — x^, T) 0 — 7r/ 2 


First, we assume a uniform wave speed 
c Q = 50 (xj — Xj) sec -1 

and compute a segment of the wave emanating from 
x 0 , 7/ 0 . Figure 3(a), in which the surface is displayed 
from the side, shows the rays and wave fronts for this 
segment. 

Next, it is assumed that the surface is no longer 
uniform, but is such that the wave speed increased 
in the x-direction according to the formula 

e(x) = c Q [ 1 + 0.008(x — x 0 )] 

Figure 3(b) shows the influence of the gradient on the 
rays that propagate from x 0 ,r) 0 initially with same 
directions as those in figure 3(a). 

Figure 4 demonstrates the effect of the local sur- 
face total curvature on ray spreading. In figure 4(a), 
the surface is the same as that used in the previ- 
ous example. It is elliptic ( K > 0) everywhere. 
Now, however, the rays do not spread from a sur- 
face source, but are initially parallel as if launched 
by grazing rays of a plane wave. Examples of surface 
waves initiated in this way are Franz and Stoneley 
waves. As the wave propagates over the surface, the 
rays converge and actually pass through a focus. 

Figure 4(b) shows a surface that is hyperbolic 
everywhere (K < 0). Its nose cross section is an 
elliptic whose major-to-minor axis ratio is 2:1, and 
the base cross section is a circle. The lofting functions 
are 


a(x) = R(x) 

0{x) = 0.375(0. liZ + 0.9 R 2 )(x f - xj 

Again, the rays are initially parallel, but the ray strip 
rapidly diverges. 

Application of Results to Wave Field 
Calculation 

As noted earlier, the phase is determined by the 
time along the ray, for harmonic motion of fixed 
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frequency. 
a by 


It is related to the ray path arc length 



The amplitude A of the surface wave is deter- 
mined according to an energy principle in reference 2 
by the formula 


Constant 

v^Go /4 


(37) 


where p is a parameter that depends on the densities 
on the two sides of the surface. If p and c are given, 
then A is fully determined once the deviation \/Ul 
of the rays has been computed by one of the two 
methods previously described. Equation (37) may 
be compared with the formula given in reference 3 
for a surface disturbance U that satisfies the wave 
equations and the impedance boundary condition 


then in the direction orthogonal to the ray (in the 
N, h) plane. 

If <f> = 7t/2 — 9 is the angle at which a launched 
ray leaves the surface at a distance a along the ray, 
then the angle at which the ray at a 4- da is launched 
is (f> + ^ da so that 


d<f> 



(41) 


If c is constant along the ray path, 9 is constant, 
and therefore results from the normal curvature 
of the ray; 

d<j> = da (42) 

However, if c varies along the ray, there is also a 
corresponding variation in 9 , which is obtained by 
differentiating equation (40) along a: 


dO _ c m dc 
da c 2 cos 9 da 


(43) 


-j- + kYU = 0 (38) 

QjTX 

where k is the reference wave number, and Y is 
proportional to the surface impedance, which may 
vary over the surface. In reference 3, Y is taken to 
be real and positive. Instead of treating a surface 
ray strip, reference 3 treats an actual tube of rays 
adjacent to the surface. The resulting amplitude 
formula is 

\/ c y 

A = Constant , (39) 

G l J 4 

If the surface wave is of a type that ripples the 
surface and if the surrounding medium will support 
a wave, then the surface wave may launch a head 
wave into this medium. If the sound speed c m in the 
medium is less than the surface wave speed, then the 
rays leave the surface at an angle 0 to the surface 
normal such that 

sin 9 = — (40) 

c 

The surface wave is “supersonic” relative to the 
surrounding medium, and the head wave is the Mach 
wave associated with that motion. If the wave speed 
c m is constant, the rays in the surrounding medium 
are straight. In this case, the variation of ray tube 
area in the medium depends only on the incremental 
differences in the ray directions as they leave the 
surface. To determine these incremental angles, we 
study a ray strip associated with a particular axial, 
or central, ray. First, we consider the variation in 
ray directions along the ray (in the N,T plane), and 


Consequently, the total incremental angle is 

d<f> = (« n l “ 2 Cm a ~f] dcr ( 44 ) 

\ c l cos 9 da J 

Next, we consider the variation in launch angle 
normal to the TV, T plane. Let p be a unit vector in 
the launch direction: 

p = sin 6T + cos ON (45) 

Denote the curvature vector of the surface wave front 
by k u . Its component in the launch direction is 

/C/ = Kfj * p 

= (/%2^V + Kg2T) ■ p 
= K n 2 cos 6 + k 9 2 sin 9 (46) 

The angular ray spreading orthogonal to the TV,T 
plane is 

dip = ki dn (47) 

Finally, the ray tube spreading for the launched 
rays is proportional to which, by using equa- 

tions (44) and (47), yields 

Spreading = Constant ( K n \ (/c n 2 cos0 + k q 2 sin0) 

\ c z cos0 dcr ) 

(48) 

In this equation, n g 2 is obtained from equation (31). 
The normal curvature K n \ along the ray direction 

e + 2/A + gA 2 

nl E + 2FX + GA 2 1 j 
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while K n 2 is computed from the same formula with X 
replaced by A n . For Franz and Stoneley waves, 6 = 
90°; therefore c = c m and ^ = 0 in equation (43), 
and equation (48) reduces to 

Spreading = Constant (ft n i^2) (50) 


Concluding Remarks 

Equations have been derived for computing sur- 
face waves, including surfaces with a nonuniform 
wave speed. The prior literature has dealt primar- 
ily with the theoretical development with little con- 
sideration given to computational methods, and ex- 
amples have been limited to waves on surfaces of 
simple analytic description such as cones, spheres, 
cylinders. The computational procedure presented 
herein is a relatively general method. Sample calcu- 
lations illustrated the procedure for a class of prac- 
tical shapes of the type that include aerodynamic 
and hydrodynamic surfaces. Equations for spreading 
of rays launched from the surface into a surround- 
ing medium that will support waves have also been 
included. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
June 9, 1989 
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Figure 2. Top view of surface used for sample calculation 


(a) Constant wave speed 


(b) Small positive wave speed gradient. 
Figure 3. Surface for sample calculation 
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